#!/usr/bin/perl -w

use strict;

use Math::Trig;
use DateTime;
use DateTime::Format::Epoch::JD;

use GD;

my $x_size = 2000;
my $y_size = 2000;

my $img = new GD::Image($x_size, $y_size);
my $black = $img->colorAllocate(0, 0, 0); 
my $white = $img->colorAllocate(255,255,255);
$img->fill(0, 0, $white);
#$img->transparent($white);
$img->interlaced('true');
$img->setAntiAliased($white);
$img->rectangle(0, 0, $x_size - 1, $y_size - 1, $black);

# coeff pour la projection sur la carte
my $coeff = 1000;

##
# Date
#

#my $date = DateTime->now;
my $date = DateTime->new(
	year => 2011, month => 8, day => 14,
	hour => 23, minute => 0,
);


##
# Coordonees Omey
##

my $latitude = "48.95N"; # Nord; positif, Sud: negatif
my $longitude = "4.5E"; # Ouest: negatif, Est: positif

my $lat = deg2rad(48.95);
my $longi = deg2rad(4.5);

##
# Calcul du jour Julien
##

my $julian = DateTime::Format::Epoch::JD->format_datetime($date);


##
# Calcul de l'heure siderale
##

#Nous connaissons le jour julien, nous en déduisons le nombre de siècle depuis le 01/01/2000 grâce à la formule suivante :
my $T = ( $julian - 2451545 ) / 36525;

#On en déduit l'heure sidérale en seconde grâce à la formule suivante :
my $H1 = 24110.54841 + ( 8640184.812866 * $T) + ( 0.093104 * ( $T^2 ) ) - (0.0000062 * ( $T^3 ) );

#en heure l'heure sidérale est donc :
my $HSH = $H1 / 3600;

#Mais il faut ramener cette heure dans un intervalle de 0 à 24H et ne garder que la partie fractionnaire de ce nombre d'où :
my $HS = (( $HSH / 24 ) - int( $HSH / 24 ))*24;

#Calcul du décalage dû à l'heure sidérale
#Sachant que la Terre tourne sur elle-même en 23H56min4s, on en déduit l'angle auquel correspond l'heure sidérale
#my $angleH = 2 * pi * $HS / (23H56min4s)
my $angleH = 2 * pi * $HS / 24;

#Calcul du décalage dû à l'heure 
#Il depend de l'heure en temps universelle, c'est à dire l'heure à Greenwich, d'où le paramètre "-zone", il est compté à partir de 12H (car l'heure sidérale est comptée à midi) d'où le paramètre "-12", ce décalage dépend aussi de la rotationde la Terre en 23H56min4s.
#angleT = (heure - 12 + minute/60 - zone) * 2 * PI/(23H56min4s)
my $angleT = ($date->hour - 12 + $date->minute/60) * 2 * pi / 24;

#print "date=$date, jour julien = $julian, T=$T, H1=$H1, HSH=$HSH, HS=$HS, angleH=$angleH, angleT=$angleT\n";


##
# Couleur des etoiles
##

my %star_colors = (
	O => $img->colorAllocate(155, 176, 255),  #9bb0ff  
	B => $img->colorAllocate(170, 191, 255),  #aabfff  
	A => $img->colorAllocate(202, 215, 255),  #cad7ff  
	F => $img->colorAllocate(248, 247, 255),  #f8f7ff  
	G => $img->colorAllocate(255, 244, 234),  #fff4ea  
	K => $img->colorAllocate(255, 210, 161),  #ffd2a1  
	M => $img->colorAllocate(255, 204, 111),  #ffcc6f  
	C => $img->colorAllocate(255, 0, 0),  # carbon stars
	NA => $img->colorAllocate(0, 0, 0),  # unknow type
);

##
# Catalogue Hipparcos
##

my %hpstars;

my $hpstars_catalog = "hp.csv";
open HPSTARS, "<$hpstars_catalog" or die "unable to open Hipparcos Catalog";
while (<HPSTARS>) {
        next if /^#/;
        chomp;
        #my ($raj2000, $dej2000, $hp, $mag, $plx, $pmRA, $pmDE, $e_Plx, $bv, @notes) = split /;/;
        my ($raj2000, $dej2000, $hp, $mag) = split /;/;
	($raj2000 = $raj2000) =~  s /\s//g;
	($dej2000 = $dej2000) =~  s /\s//g;
	($hp = $hp) =~  s /\s//g;
        $hpstars{$hp}->{ra} = $raj2000; # RA en H:M:S
        $hpstars{$hp}->{dec} = $dej2000; # DEC en degres
        $hpstars{$hp}->{mag} = $mag;
}
close HPSTARS;

##
# Catalogue Brightstars
##

my %stars;

my $stars_catalog = "brightstars.csv";
open STARS, "<$stars_catalog" or die "unable to open Brighstars Catalog";
while (<STARS>) {
        next if /^#/;
        chomp;
        my ($raj2000, $dej2000, $hr, $name, $hd, $ads, $varid, $ra2000, $de2000, $mag, $bv, $sptype, @notes) = split /;/;
	($hd = $hd) =~  s /\s//g;
        $stars{$hd}->{ra} = $raj2000; # RA en H:M:S
        $stars{$hd}->{dec} = $dej2000; # DEC en degres
        ($stars{$hd}->{name} = $name) =~ s/\s//g;
        $stars{$hd}->{mag} = $mag;
	my $color = $1 if $sptype =~ /([OBAFGKM])/;
        $stars{$hd}->{color} = $color || "NA";
        $stars{$hd}->{desc} = join ', ', @notes;
}
close STARS;

for my $i (sort { $stars{$a} <=> $stars{$b} } keys %stars) {
	next if $stars{$i}->{mag} > 6;

	my ($hau, $X, $Y) = coords($stars{$i}->{ra}, $stars{$i}->{dec});
	
	next if $hau < 0;

	my $mag = $stars{$i}->{mag};
	my $size = 10 - $mag;
	$img->arc($X, $Y, $size, $size, 0, 360, $star_colors{$stars{$i}->{color}});
	#$img->arc($X, $Y, $size, $size, 0, 360, $black);
	$img->fill($X, $Y, $star_colors{$stars{$i}->{color}});

	my $name = $stars{$i}->{name};
	$img->string(gdSmallFont, $X + 10, $Y, $name, $black) if $mag < 3;
	#print "$i: name=$name, mag=$mag\n" if $mag < 3;
}

##
# Catalogue Carbon Stars
##

my %cstars;

my $cstars_catalog = "carbon.csv";
open CSTARS, "<$cstars_catalog" or die "unable to open Carbon Stars Catalog";
while (<CSTARS>) {
        next if /^#/;
        chomp;
        my ($raj2000, $dej2000, $number, $mag) = split /;/;
        $cstars{$number}->{ra} = $raj2000; # RA en H:M:S
        $cstars{$number}->{dec} = $dej2000; # DEC en degres
        $cstars{$number}->{mag} = $mag;
        $cstars{$number}->{color} = "C";
}
close CSTARS;

for my $i (sort { $cstars{$a} <=> $cstars{$b} } keys %cstars) {
	next if $cstars{$i}->{mag} > 7;

	my ($hau, $X, $Y) = coords($cstars{$i}->{ra}, $cstars{$i}->{dec});
	
	next if $hau < 0;

	my $mag = $cstars{$i}->{mag};
	my $size = 10 - $mag;
	$img->arc($X, $Y, $size, $size, 0, 360, $star_colors{$cstars{$i}->{color}});
	$img->fill($X, $Y, $star_colors{$cstars{$i}->{color}});
}


##
# Constellations lines
##

my $lines_data = "constellationship.fab";
open LINES, "<$lines_data" or die "unable to open constellations lines data";
while (<LINES>) {
	chomp;
	my ($const, $nb, @points) = split /\s+/;

	#print "const=$const, nb=$nb, nbpoints=$#points, points=", join(':', @points), "\n";

	while ($#points > 0) {
		my $star1 = pop @points;
		my $star2 = pop @points;

		my ($s1_hau, $s1_X, $s1_Y) = coords($hpstars{$star1}->{ra}, $hpstars{$star1}->{dec});
		my ($s2_hau, $s2_X, $s2_Y) = coords($hpstars{$star2}->{ra}, $hpstars{$star2}->{dec});

		#print "\tstar1=$star1 (", $hpstars{$star1}->{ra}, ",", $hpstars{$star1}->{dec}, "): hau=$s1_hau, X=$s1_X, Y=$s1_Y\n";
		#print "\tstar2=$star2 (", $hpstars{$star2}->{ra}, ",", $hpstars{$star2}->{dec}, "): hau=$s2_hau, X=$s2_X, Y=$s2_Y\n";

		if ($s1_hau > 0 && $s2_hau > 0) {
			$img->line($s1_X, $s1_Y, $s2_X, $s2_Y, $black)
		}
	}
}

##
# Constellations boundaries
#

#Column _RAJ2000        (F9.5)  Right ascension (FK5) Equinox=J2000. (computed by VizieR, not part of the original data)        [ucd=pos.eq.ra;meta.main]
#Column _DEJ2000        (F9.5)  Declination (FK5) Equinox=J2000. (computed by VizieR, not part of the original data)    [ucd=pos.eq.dec;meta.main]
#Column RAdeg   (F10.6) Right ascension in decimal hours (J2000)        [ucd=pos.eq.ra;meta.main]
#Column DEJ2000 (F11.7) Declination in degrees (J2000)  [ucd=pos.eq.dec;meta.main]
#Column cst     (A4)    Constellation abbreviation      [ucd=meta.id.part]
#Column type    (A1)    [OI] Type of point (Original or Interpolated)   [ucd=meta.note]


#my $constb_catalog = "constellations_boundaries.csv";
#open CONSTB, "<$constb_catalog" or die "unable to open constellations boundaries";
#while (<CONSTB>) {
#        chomp;
#        my ($raj2000, $dej2000, $ra2000, $de2000, $cst, $type) = split /;/;
#        next if /^#/;
#
#	next if $cst ne "AND ";
#
#	my ($hau, $X, $Y) = coords($stars{$i}->{ra}, $stars{$i}->{dec});
#
#	next if $hau < 0;
#}
#
#close CONSTB;


##
# Catalogue NGC
##

my %ngc;

my $ngc_catalog = "ngc.csv";
open NGC, "<$ngc_catalog" or die "unable to open NGC Catalog";
while (<NGC>) {
	chomp;
	my ($raj2000, $dej2000, $name, $type, $ra2000, $de2000, $source, $const, $size, $mag, @desc) = split /;/;
	next if /^#/;
	next if $name =~ /^I/;
	($name = $name) =~  s /\s//g;
	#($raj2000 = $raj2000) =~  s /\s//g;
	$ngc{$name}->{ra} = $raj2000; # RA en H:M:S
	$ngc{$name}->{dec} = $dej2000; # DEC en degres
	$ngc{$name}->{type} = $type;
	$ngc{$name}->{const} = $const;
	$ngc{$name}->{size} = $size;
	$ngc{$name}->{mag} = $mag;
	$ngc{$name}->{desc} = join ', ', @desc;
}
close NGC;

# sort by RA
for my $i (sort { $ngc{$a}->{ra} <=> $ngc{$b}->{ra} } keys %ngc) {
	next unless $ngc{$i}->{desc} =~ / = M/;

	#next unless $ngc{$i}->{type} =~ /Pl/;

	my$name = "NGC$i";
	$name = $1 if $ngc{$i}->{desc} =~ / = (M\d+)/;

	my ($hau, $X, $Y) = coords($ngc{$i}->{ra}, $ngc{$i}->{dec});

	next if $hau < 0;

	my $type = $ngc{$i}->{type};

	$img->string(gdMediumBoldFont, $X, $Y, "+" . $type . " " . $name, $black);

	print "NGC$i name=$name, type=$type, X=$X, Y=$Y, hau=$hau, RA=", $ngc{$i}->{ra}, ", DEC=", $ngc{$i}->{dec}, ", desc=", $ngc{$i}->{desc}, "\n";
}

$img->arc($x_size / 2, $y_size / 2 , $x_size - 20 , $y_size - 20 , 0, 360, $black);

open IMG, ">chart.png";
print IMG $img->png;
close IMG;


sub coords() {
	my $asc = deg2rad(shift);
	my $dec = deg2rad(shift);

	my $angle = $angleH + $angleT;

	my $H = $angle - $asc + $longi;

	##
	# Calcul de la hauteur
	##

	my $sinushauteur = sin($dec) * sin($lat) - cos($dec) * cos($lat) * cos($H);

	my $hau = asin ($sinushauteur);

	##
	# Calcul de l'azimut
	##

	my $cosazimuth = ( sin($dec) - sin($lat) * sin($hau) ) / ( cos($lat) * cos($hau) );
	#L'azimut est un angle compris entre 0 et 360°, nous avons donc besoin d'un calcul intermédiaire :
	my $sinazimuth = ( cos($dec) * sin ($H) ) / cos($hau);
	my $az;
	if ($sinazimuth > 0) {	
		$az = acos($cosazimuth);
	}
	else {
		$az = - acos($cosazimuth);
	}

	my $dist = -$coeff * ( (-2/pi) * $hau + 1 );
	#(le paramètre - permet de placer le Nord en haut de la carte, le coeff permet de déterminer la taille de la carte et peut varier a volonter) Sur les axes X et Y, on a donc :
	my $X = $dist * cos($az) + $x_size / 2;
	my $Y = $dist * sin($az) + $y_size / 2;
	
	return ($hau, $X, $Y);
}

__END__

=pod

=head1 Les formules

voir http://emilie.bodin.free.fr/logiciel/logiciel.html

=head2 Notations utilisées

Latitude = lat
Longitude = longi
Déclinaison = dec
Ascension droite = asc
Azimut =az
Hauteur = hau
Angle Horaire de l'étoile = H = angle - asc + longi
angleH = angle lié a l'heure sidérale.
angleT = angle lié a l'heure.
angle = angleH + angleT

=head2 Calcul de la Hauteur

sinushauteur = sin(dec) * sin(lat) - cos(dec) * cos(lat) * cos(H)
La hauteur est un angle compris entre -90° et +90°, la hauteur s'obtient donc simplement par :
hau = arcsin (sinushauteur)

Calcul de l'Azimut :

cosazimuth = ( sin(dec) - sin(lat) * sin(hau) ) / ( cos(lat) * cos(hau) )
L'azimut est un angle compris entre 0 et 360°, nous avons donc besoin d'un calcul intermédiaire :
sinazimuth = ( cos(dec) * sin (H) ) / cos(hau)
Si sinazimuth > 0 alors :
az = + arccos(cosazimuth)
Sinon :
az = - arccos(cosazimuth)

=head2 Calcul du jour julien:

Pour la date, il faut utiliser un système de référence plus simple que le système jour-mois-année traditionnel, on utilise en général le JOUR JULIEN : ce calendrier est très simple, il ne fait que compter les jours à partir d'une date de référence. La date de référence est le 1er janvier de l'an -4712 à 12H00 (par exemple le 1er janvier 2OOO à 00H00 correspond au jour julien 2451544.5). 
L'utilisateur a donc entré une Année, un Mois, un Jour, une Heure et une Zone, nous en déduisons par le calcul le nombre de jour depuis le 1er janvier 2000 auquel il correspond. L'heure donnée est l'heure locale, la zone correspond au décalage par rapport à Greenwich.

[ Par convention : la fonction Int est la fonction qui rend la partie entière d'un nombre ]

Si le mois est inférieur à 3 :
  alors Mois = Mois + 12 et Année = Année - 1
  sinon on ne fait rien


Puis on fait les calculs suivants :
A = Int( Année / 100 )
B = 2 - A + Int( A / 4 )

( les termes A et B sont des termes correctifs qui doivent être utilisés pour les dates à partir du 15 octobre 1582, date de la reforme du calendrier julien en calendrier grégorien. Dans mon logiciel, j'ai supposé que l'utilisateur tape toujours une date apres celle-ci )
C = Int( 365.25 * Année )
D = Int( 30.6001 * ( Mois + 1 ) )
Jour Julien :JJ = B + C + D + jour + 1720994.5


=head2 Calcul de l'heure sidérale

L'heure sidérale va nous permettre de calculer un angle qui dépend de l'heure d'observation et de la date.
Nous connaissons le jour julien, nous en déduisons le nombre de siècle depuis le 01/01/2000 grâce à la formule suivante :
T = ( JJ - 2451545 ) / 36525
On en déduit l'heure sidérale en seconde grâce à la formule suivante :
H1 = 24110.54841 + ( 8640184.812866 * T) + ( 0.093104 * ( T^2 ) ) - (0.0000062 * ( T^3 ) )
en heure l'heure sidérale est donc :
HSH = H1 / 3600
Mais il faut ramener cette heure dans un intervalle de 0 à 24H et ne garder que la partie fractionnaire de ce nombre d'où :
HS = (( HSH / 24 ) - Int( HSH / 24 ))*24
Calcul du décalage dû à l'heure sidérale

Sachant que la Terre tourne sur elle-même en 23H56min4s, on en déduit l'angle auquel correspond l'heure sidérale
angleH = 2 * PI * HS / (23H56min4s)

Calcul du décalage dû à l'heure 

Il depend de l'heure en temps universelle, c'est à dire l'heure à Greenwich, d'où le paramètre "-zone", il est compté à partir de 12H (car l'heure sidérale est comptée à midi) d'où le paramètre "-12", ce décalage dépend aussi de la rotationde la Terre en 23H56min4s.
angleT = (heure - 12 + minute/60 - zone) * 2 * PI/(23H56min4s)


Projection sur un plan :

Pour obtenir une représentation des étoiles, il ne reste plus qu'à projeter ces coordonnées sur un plan. 
Il ne faut pas oublier que seul les étoiles donc la hauteur (hau) est positive sont visibles ! (les autres sont au-dessous de nos pieds).
Au départ, on pourrait penser faire une projection de la carte avec un cos(hau) comme on le ferait classiquement pour projeter une sphère sur un plan en mathématique, en fait une projection linéaire donne un meilleur résultat et une carte moins déformée sur les bords. Nous utilisons donc :
dist = -coeff * ( (-2/PI) * hau + 1 )
(le paramètre - permet de placer le Nord en haut de la carte, le coeff permet de déterminer la taille de la carte et peut varier a volonter) Sur les axes X et Y, on a donc :
X = dist * cos(az)
Y = dist * sin(az)

Autres types de projection :

http://sam.electroastro.pagesperso-orange.fr/dossiers/projection/project.htm


Objets visibles par saison (tries par RA):
0 to 3 => Oct
3 to 6 => Nov/Dec
6 to 9 => Jan
9 to 12 => Feb/Mar
12 to 15 => Apr
15 to 18 => May/June
18 to 21 => July
21 to 0 => Aug/Sept

=cut
